function [ out ] = fun_d2sigma_new2(theta,argument,m_regr,combos_all_unique,J,ind1,ind2)

B_der = zeros(1,size(combos_all_unique,1));

for i=1:size(combos_all_unique,1)
    
    ind_other = perms(1:J-1);
    
    combos_x_other = combos_all_unique(i,2:J);
    combos_x_other_r = combos_x_other(ind_other);
    
    combos_z_other = combos_all_unique(i,J+2:2*J);
    combos_z_other_r = combos_z_other(ind_other);
    
    v = ([repmat(combos_all_unique(i,1),[size(ind_other,1),1])   combos_x_other_r...
        repmat(combos_all_unique(i,J+1),[size(ind_other,1),1]) combos_z_other_r]) - 1;
    
    v = unique(v,'rows');       % KEY to avoid counting multiple times the combos that correspond to only one of combos_all
    
    v_minus1_minus1 = v;
    v_minus1_minus1(:,[ind1 ind2]) = v_minus1_minus1(:,[ind1 ind2])-1;
    
    v_plus1_plus1 = v;
    v_plus1_plus1(:,[ind1 ind2]) = v_plus1_plus1(:,[ind1 ind2])+1;
    
    v_minus1_0 = v;
    v_minus1_0(:,ind1) = v_minus1_0(:,ind1)-1;
    
    v_plus1_0 = v;
    v_plus1_0(:,ind1) = v_plus1_0(:,ind1)+1;
    
    v_0_minus1 = v;
    v_0_minus1(:,ind2) = v_0_minus1(:,ind2)-1;
    
    v_0_plus1 = v;
    v_0_plus1(:,ind2) = v_0_plus1(:,ind2)+1;
    
    v_minus1_plus1 = v;
    v_minus1_plus1(:,ind1) = v_minus1_plus1(:,ind1)-1;
    v_minus1_plus1(:,ind2) = v_minus1_plus1(:,ind2)+1;
    
    v_plus1_minus1 = v;
    v_plus1_minus1(:,ind1) = v_plus1_minus1(:,ind1)+1;
    v_plus1_minus1(:,ind2) = v_plus1_minus1(:,ind2)-1;
    
    v_0_0 = v;
    
    B_der(i) = sum((prod(BernsteinElement_2(v_minus1_minus1,m_regr,argument'),3).*(m_regr(ind1)-(v(:,ind1)-1)).*(m_regr(ind2)-(v(:,ind2)-1)) +...
        prod(BernsteinElement_2(v_0_minus1,m_regr,argument'),3).*(2*v(:,ind1)-m_regr(ind1)).*(m_regr(ind2)-(v(:,ind2)-1)) +...
        -prod(BernsteinElement_2(v_plus1_minus1,m_regr,argument'),3).*(v(:,ind1)+1).*(m_regr(ind2)-(v(:,ind2)-1)) +...
        prod(BernsteinElement_2(v_minus1_0,m_regr,argument'),3).*(m_regr(ind1)-(v(:,ind1)-1)).*(2*v(:,ind2)-m_regr(ind2)) +...
        prod(BernsteinElement_2(v_0_0,m_regr,argument'),3).*(2*v(:,ind1)-m_regr(ind1)).*(2*v(:,ind2)-m_regr(ind2)) +...
        -prod(BernsteinElement_2(v_plus1_0,m_regr,argument'),3).*(v(:,ind1)+1).*(2*v(:,ind2)-m_regr(ind2)) +...
        -prod(BernsteinElement_2(v_minus1_plus1,m_regr,argument'),3).*(m_regr(ind1)-(v(:,ind1)-1)).*(v(:,ind2)+1) +...
        -prod(BernsteinElement_2(v_0_plus1,m_regr,argument'),3).*(2*v(:,ind1)-m_regr(ind1)).*(v(:,ind2)+1) +...
        +prod(BernsteinElement_2(v_plus1_plus1,m_regr,argument'),3).*(v(:,ind1)+1).*(v(:,ind2)+1))',2);
    
end

out = B_der*theta;

end